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Recent theories for supersonic jet noise have used an instability wave noise generation 
model to predict radiated noise. This model requires a known mean flow that has typically 
been described by simple analytic functions for single jet mean flows. The mean flow of 
supersonic coaxial jets is not described easily in terms of analytic functions. To provide 
these profiles at all axial locations, a numerical scheme is developed to calculate the mean 
flow properties of a coaxial jet. The Reynolds-averaged, compressible, parabolic boundary 
layer equations are solved using a mixing length turbulence model. Empirical correlations 
are developed to account for the effects of velocity and temperature ratios and Mach 
number on the shear layer spreading. Both normal velocity profile and inverted velocity 
profile coaxial jets are considered. The mixing length model is modified in each case to 
obtain reasonable results when the two stream jet merges into a single fully developed jet. 
The mean flow calculations show both good qualitative and quantitative agreement with 
measurements in single and coaxial jet flows. 
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1. INTRODUCTION 

The current interest in the development of a High Speed Civil Transport that is 
economically viable and environmentally compatible has renewed research efforts on 
supersonic jet mixing and noise reduction. Supersonic jets are intense noise generators, and 
methods must be found to modify the noise generation process to reduce radiated noise 
levels if a future high speed civilian aircraft is to meet community noise regulations. In 
a recent review, Seiner and Krejsa [1] discuss the challenge of reducing supersonic jet noise 
associated with both mixing and shocks while maintaining acceptable propulsion system 
performance. The goal is to maintain the total thrust and mass flow of the propulsion 
system at the highest levels while reducing the strength of the noise sources. 

One method to modify the noise generated by a supersonic jet is to replace the single 
stream jet with a dual stream, coaxial jet. Each stream may have a different initial velocity 
and possibly different initial temperatures. For classification purposes, when the coaxial 
jet flow has a higher inner stream velocity than an outer stream velocity, it is referred to 
as a normal velocity profile (NVP) jet, as illustrated in Figure 1(a). If the outer stream 
velocity is higher than the inner stream velocity, the jet has an inverted velocity profile 
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(I VP); see Figure 1(b). The temperature of the streams is arbitrary but, for practical 
propulsion systems, the higher velocity stream typically has the higher temperature. 

The analysis to predict the noise generated from supersonic, perfectly expanded jets is 
based on a theory that instability waves or large scale coherent structures control the 
development of free jet shear flows. When the instability wave phase velocities exceed the 
ambient speed of sound, these waves are the dominant source of mixing noise radiated 
into the downstream arc of the jet. Even when shocks are present in the jet flow, the noise 
from the jet that radiates into the downstream arc is due primarily to turbulent mixing, 
whereas the broadband shock associated noise dominates in the upstream arc. Extensive 
measurements have shown this for single jets [2] and for coaxial jets [3]. The assumption 
that the jet is perfectly expanded simplifies the analysis and allows us to concentrate on 
profile shaping as a means to further reduce the mixing noise. 

To complete the analysis in the instability wave noise generation model, the equations 
of motion are decomposed with the assumption that the flow variables are a combination 
of a mean flow component and a fluctuating component representing the large scale 
coherent structures. This results in equations that describe the instability waves or large 
scale structures that are superposed on a known mean flow field. Thus, the mean flow is 
an input to the instability wave problem and must be determined separately. A complete 
solution to the instability noise generation model requires that the developing mean flow 
properties be known at every axial location. In the present paper, we describe an analysis 
to determine the developing mean flow properties of a compressible coaxial jet. Due to 
the complex nature of developing NVP and IVP jets with various velocity and temperature 
conditions, the mean flow has been determined numerically. In related papers, we will 



X 

Figure I. Illustrations of coaxial jet flows: (a) normal velocity profile coaxial jet; (b) inverted velocity profile 
coaxial jet. 
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describe the evolution of the instability waves or large scale structures and the noise that 
they generate and radiate [4, 5]. 

For the most part, in the past, in calculations of instability waves in a single 
axisymmetric jet, analytic functions have been used to characterize the mean flow. These 
analytic functions have been based on results from experimental measurements where data 
have been correlated by local similarity variables. Examples of this approach include the 
measurements of: Lau et al. [6j, whose data were fitted by an error function profile; Cohen 
and Wygnanski [7], who fitted data with a series of hyperbolic tangent functions, to add 
corrections to the classical hyperbolic tangent profile; and Morrison and McLaughlin [8], 
who fitted data by a haif-Gaussian profile. Michalke [9] has summarized the use of different 
analytic functions in the calculation of instability waves. The measured data typically 
include only velocity profiles, which are sufficient for incompressible instability wave 
calculations. When compressibility or mean density effects are important, the instability 
wave calculations require that either the temperature or the density profile be specified. 
Often, the approximations have been made that allow the Crocco-Busemann relation to 
be used. This describes the temperature or density profile as a function of the velocity alone 
[9-1 1]. Using the methods to be described in this paper, Dahl [12] calculated velocity and 
density profiles for both subsonic and supersonic jets and compared the results to 
hyperbolic tangent velocity profiles and the corresponding density profiles from the 
Crocco-Busemann relation. Good agreement was found for both isothermal and hot jets 
showing the usefulness of the Crocco-Busemann relation in single jet instability wave 
calculations. With the mean profile information defined analytically, the mean profile 
variables and all their derivatives are available at all points in the radial direction for the 
numerical calculation of instability wave characteristics. Furthermore, these analytic 
functions may be continued into the complex plane. This allows for the calculation of 
damped in viscid waves [13]. Thus, the ability to use analytic functions to describe the mean 
flow simplifies the study of single jet instability waves. 

In the case of a coaxial jet, it is not clear how to develop analytic functions to describe 
the profiles at all axial locations. Tam and Burton [14], in a single supersonic jet case, used 
a generalized half-Gaussian function to describe the mean velocity at all axial locations. 
The scaling parameters, the centerline velocity, the core radius, and the half-width of the 
mixing region, were defined by cubic spline fits to measured data. The density profile was 
found by keeping the total temperature constant. This approach was possible due to the 
availability of measured data; however, there is little measured data for coaxial jets, 
especially with supersonic conditions, that would allow an analytical description to be 
made at all axial locations, including the merging region of a normal profile and an 
inverted profile into a single jet. Furthermore, whereas temperature profiles have been 
measured in a subsonic jet, they are not typically measured in a supersonic jet. Therefore, 
even though velocity measurements may be available for coaxial jets with supersonic 
conditions, it is doubtful that the Crocco-Busemann relationship could provide the 
appropriate temperature profiles for either normal or inverted conditions. Thus, in the 
present paper the mean profiles for coaxial jets are determined numerically. Morris [15, 16] 
has calculated instability waves using numerically generated velocity profiles for a single 
incompressible jet. The significant extension here is the inclusion of compressibility effects 
and the consideration of both normal and inverted velocity profile coaxial jets. 

Techniques to calculate the mean flow properties of coaxial jets numerically have ranged 
from the use of Reichardt’s theory for free turbulence mixing [17-19] to the use of full 
Navier-Stokes solvers [20]. Early studies primarily modelled subsonic jets and used simple 
formulations in the turbulence models for the empirically derived constants. More recent 
studies have used two equation k—E turbulence models. In all cases, the calculated velocity 
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profiles resembled measured data with varying degrees of accuracy. However, many of 
these numerical techniques lack an ease and robustness that is necessary to predict the 
development of a variety of NVP and l VP jets with a wide variety of operating conditions. 
To obtain these qualities, a numerical technique based on the approach of Patankar and 
Spalding [21] has been used. In their work, the development of a turbulent, axisymmetric 
free jet is assumed to be governed by the Reynolds-averaged, compressible, parabolic 
boundary layer equations. A turbulence model is required to complete the set of equations. 
In the next section, we describe the governing boundary layer equations and discuss the 
basic turbulence model for a single compressible shear layer with different velocity and 
density ratios. The modifications to the turbulence model are then described for the 
merging flow in NVP and IVP jets. Next, the numerical formulation is presented. Finally, 
numerical predictions are made for both single and coaxial jets that compare calculated 
profiles to measured data and other calculated results for the same jet operating conditions. 
It should be noted that the goal of this part of the work is the calculation of the mean 
flow properties of a coaxial jet, either with a normal initial velocity profile or with an 
inverted initial velocity profile, that expands and merges in a smooth and physically 
reasonable fashion in the downstream direction. It should be recalled that the resulting 
predictions of the mean flow development form the basis for an instability wave model 
of the large scale coherent structures and their noise radiation. Therefore, the absolute 
accuracy of the calculated mean flows compared to any particular set of measured data 
is secondary to whether or not the calculated results bear a strong qualitative resemblance 
to the behavior of coaxial jet flows. However, it will be shown that the predictions, 
when compared to available measurements, actually show reasonable quantitative 
agreement. 


2. MEAN FLOW DEVELOPMENT 
2.1. EQUATIONS FOR COMPRESSIBLE flow 

The compressible flow of a free jet is governed by the continuity, momentum, energy 
and state equations. Following the discussion in Anderson et al. [22], we assume that the 
mean flow development of an axisymmetric free jet is governed by the Reynolds-averaged 
boundary layer equations for compressible flow; since the flow is also assumed to be 
turbulent. The flow variables are decomposed into a time averaged or mean value part and 
a fluctuating or turbulent part (e.g., p = p + p\ where the overbar denotes the time 
averaged part and the prime denotes the fluctuating part). With the jet static pressure 
matched to the ambient pressure and the density-velocity and density-enthalpy 
correlations neglected [23], the equations in axisymmetric co-ordinates reduce to the 
following: 

continuity, 

(rpu) 4- Y r ( r &) = °> 0) 

momentum, 





da 

^Tr~ puv )\ 


(2) 
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energy. 


where v is a mass averaged quantity, defined as 

v = (l/p)(pv 4- p'v'). (4) 

These equations have been non-dimensionalized by the following reference values: spatial 
co-ordinates by R\, velocity by Ui, density by p it pressure by p\U], enthalpy and c p T by 
U\, and viscosity by R l p ] U l . The subscript 1 is used to indicate jet exit conditions for a 
single jet or for the inner stream of a coaxial jet. 

We assume a perfect gas. Thus, the equation of state is 

p = pR g f , (5) 

where R g is the gas constant. Also, the laminar viscosity is related to the mean temperature 
according to Sutherland’s Law: 
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5*0 is 1 10-3 and po is the viscosity at the reference temperature T 0 of 273 K. 


2.2. TURBULENCE MODEL 

In equa tions (2) and (3), the Reynolds stress term —pu'v' and the Reynolds heat flux 
term —pc p v'T' are undefined. These terms must be modelled in order to solve the system 
of equations. In anticipation of the difficulties in the definition of a turbulence model that 
would work for a supersonic coaxial jet with both normal or inverted initial velocity 
profiles and streams of different temperature (an example of the difficulties in turbulence 
modelling for single stream supersonic jets is found in Dash [24]) it was decided to use 
the most basic model; namely, the mixing length model for the turbulence closure. Any 
turbulence model chosen for these calculations would require modification and calibration 
to obtain reasonable results for both normal and inverted profile coaxial jets. The mixing 
length model represents the turbulent stress from which an effective viscosity is defined. 
Through the use of a turbulent Prandtl number to relate eddy diffusivities for heat and 
momentum, an effective Prandtl number may be calculated from the effective viscosity. 

In the mixing length model, the Reynolds stress term is written as [10, 22] 

- / 7 00 /n\ 

pu v = pij- r ( 7 ) 


and the heat flux term as 


- pcpv'T ' ~ 


CpPvOT 

Pr T Or’ 


( 8 ) 


The turbulent Prandtl number Pr T has essentially a constant value for free jet calculations 
[10, 22, 25, 26]. The mixing length model is embodied in the turbulent viscosity, 

p T = p(C l C : l) 2 \du/dr\, (9) 

where l is a characteristic length scale, and Ci and C : are mixing length constants; however, 
for this model, they are taken to be functions of the flow conditions. 
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It is known that the growth of a shear layer depends on the How conditions at the shear 
layer edges with the velocities and densities of the two streams being of primary interest 
[10, 25, 27, 28]. This dependence occurs whether the flow is compressible or incompressible. 
In an attempt to separate out purely compressible effects, Papamoschou and Roshko [28] 
and others have normalized the compressible shear layer growth rate by an incompressible 
growth rate value at the same edge conditions. This followed the earlier work of Brown 
and Roshko [27], that showed the effects of different velocity and density edge conditions 
on shear layer growth for incompressible flows. Thus, in order to incorporate these effects 
into the mixing length model, to obtain the proper shear layer growth, two factors have 
been developed that depend on the flow conditions; one for incompressible conditions and 
the other for compressibility effects. 

The factor Ci is the incompressible part of the mixing length constant. It depends on 
the velocity ratio lhjU\ and the density ratio p 2 /pi between the two streams on either side 
of the shear layer, where U\ and p x are the mean velocity and density of the inner stream 
and lh and p 2 are the mean velocity and density of the outer stream. Correlations for the 
expected vorticity thickness growth rate have been developed by many investigators from 
experimental evidence. Thus, given that we know the vorticity thickness growth rate for 
a fixed velocity ratio and a fixed density ratio, we have adjusted C i (C 2 = l) until the 
calculated initial vorticity thickness growth rate of a single jet agrees with the correlated 
value. Continuing this process for a range of velocity and density ratio values results in 
a series of calibration curves for C\. 

C 2 is the compressible part of the mixing length definition. Its purpose is to decrease 
the growth of the shear layer as compressibility effects become important. It depends on 
a Mach number in a relative frame of reference convecting with the real phase speed of 
a growing disturbance in the shear layer. This convective Mach number depends on the 
velocity ratio, on the density ratio, and on the Mach number of one of the streams. The 
convective Mach number has been used to correlate normalized measured growth rates 
of compressible shear layers. We have developed an equation that fits this data in order 
to predict shear layer growth rates for given flow conditions. Thus, we have proceeded to 
determine C 2 given that Ci is allowed to take on its previously calibrated value for the given 
velocity and density ratios. C 2 is adjusted until the calculated initial vorticity thickness 
growth rate agrees with the predicted value. For this case, a single calibration curve has 
been generated. Further details on the evaluation of C\ and C 2 are given in Appendix A. 

There are several choices that may be made for the characteristic length scale for a jet 
shear layer. Since the flow is turbulent and the motion is dominated by large structures 
or instability waves, an appropriate characteristic length scale would be [27, 29] 

/ = Au/\du/dr\ max , (10) 

where M is the velocity difference across the shear layer. Defined in this way, / can be 
interpreted as a vorticity thickness. If necessary, it can be related to other length scales 
if an appropriate functional shape for the velocity profile is chosen. For example, 
/ = 0-978/w where / w is defined as the shear layer width from the location at which 
(m — m 2 )/zJu = 0-9 to the location at which ( ii — u 2 )/Aii = 01; n is described by an error 
function, and w 2 is the normalized outer stream velocity. 

2.3. COAXIAL JET MIXING LENGTH MODEL 

Initially, a coaxial jet, either with an initial normal velocity profile or an inverted velocity 
profile, has two distinct shear layers with uniform flow conditions at both edges of both 
shear layers. As a result, the mixing length model, equation (9), gives separate constant 
values Ci, C 2 , and / for each layer. The values of C and C 2 depend on the velocity ratio. 
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Figure 2. Definitions of the merging norma! velocity profile mixing length model factors: (a) initial merging 
profile; (b) advanced merging profile; (c) fully merged profile. 


the density ratio and the convective Mach number determined from the separate edge 
conditions. The characteristic length / is determined for each shear layer from equation 
(10). The turbulent viscosity p T is then calculated across the coaxial jet using local values 
of p and \du/dr\. The separate shear layers have different constants (CY, C 2 , l) so that, at 
some point between the two shear layers, the pr profile must switch from one set of 
constants to the other. As a result, an abrupt change in the p T profile occurs that is not 
desirable numerically. Initially, \dufdr\ is near zero between the two shear layers and the 
effect of the abrupt change is negligible. As the calculations progress downstream, \dujdr\ 
between the two shear layers increases and the jump in p T is appreciable. To overcome 
this problem, a hyperbolic tangent function is used to transpose from one set of constants 
to the other. This works until the outer jet core is about to disappear and the shear layers 
are about to merge. At that point, the mixing length model must be altered and the normal 
and the inverted velocity profile cases are treated separately. 

2.3.1. Normal velocity profile mixing length model 
The normal velocity profile mixing length model has been developed by observations 
of the behavior of the w and dii/dr profiles as the shear layers merge and by a comparison 
of calculations to measured data taken from Lau [30]. The merging shear layer profile 
contains a local minimum in \dujdr\ that disappears as the shear layer merges fully into 
a single jet; see Figure 2. This local minimum is used as a point to separate the two merging 
shear layers. A single characteristic length for the merging shear layers is defined as 

/ = Au max /\dujdr\ 01) 

where Au max is the larger Aii between the two values determined by using the separation 
point. The values of C, and C 2 are also determined from the edge conditions that give 
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Au„„ lx . The maximum gradient |rt7/frU„ is the largest value of Idw/fVI that occurs in the 
merging profile. This method for the determination of G. G and / for a merging normal 
profile has the advantage that, as the flow transposes into a single jet profile, the G. G, 
and / values transition into the appropriate form for a single jet. 


2.3.2. Inverted velocity profile mixing length model 

When an inverted profile starts to merge, a local maximum occurs in the velocity profile 
yielding du/dr = 0; see Figure 3(a). This point is used to identify the separation point 
between the two shear layers. As long as the inner core exists, the two merging shear layers 
are treated separately but their constants are added as follows: 

(GG /)((>/<// — (G Cll) inner + (G C 2.1) outer- 02) 

This increases /r T across the profile to mimic the increased turbulent action as the inverted 
profile starts to merge. When the inner core ends, equation (12) is no longer used and it 
is assumed that the mixing process in the outer shear layer dominates the flow, so that, 
(C\Cil), 0 ,at = (GG/W; see Figure 3(b). As with the normal profile, this later usage of G, 
G, and / transitions into the proper usage for a single jet downstream. It should be noted 
that the mixing length model gives p T = 0 at the local maximum, which is unrealistic. The 
simple solution taken here is to smooth the \dujdr\ profile, and hence smooth p r [12]. 


U mnr 





Figure 3. Definitions of the merging inverted velocity profile mixing length model factors: (a) initial merging 
profile; (b) advanced merging profile, no inner stream potential core; (c) fully merged profile. 
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2.4. NUMERICAL FORMULATION 

Many numerical techniques are available for the solution of the boundary layer 
equations. Any of these could have provided adequate mean flow results; but our criterion 
for mean flow calculations dictated that the mean flow results be sufficiently resolved in 
order to be used for instability wave calculations. Furthermore, we wished to avoid any 
numerical calculation of v. Hence, initially, we followed the stream function approach used 
by Patankar and Spalding [21] and Crawford and Kays [26]. 

The equations of motion are transformed to stream function co-ordinates using 

rpw = dlP/dr, rpv — —dW /dx. (13) 

These equations ensure that the continuity equation (1) is satisfied. After substitution of 
equations (7) and (8) into equations (2) and (3) and with the use of the definition for total 
enthalpy H — c p T + u 2 j2, followed by a transformation from (x, r) co-ordinates to (x, IP) 
co-ordinates, we obtain for momentum and energy: 



where 


Peff ~ P-V PJ, 


Piv = 


1 + t*rln 
1/Pr + p T /Pr T p 


(14) 

(15) 


(16, 17) 


These equations can be differenced in a variety of ways: explicitly [31], implicitly on an 
evenly spaced IP -grid [32], or transformed to a normalized IP-grid and implicitly 
differenced as in Patankar and Spalding [21]. Each of these numerical methods have been 
found to have problems. The explicit DuFort-Frankel method [31] has stability problems; 
the implicit Crank-Nicolson method using an evenly spaced IP -grid [32] does not provide 
sufficient resolution at the outer, low speed edges of the jet; and Patankar and Spalding’s 
method [21] has difficulties with entrainment boundary conditions at the outer edge. Each 
of these problems have been overcome by using fully implicit differencing and what is 
considered to be a natural grid stretching and a natural outer boundary entrainment. By 
choosing a fully implicit method, the numerical problem is inherently stable. The problem 
of grid resolution is overcome by the use of an evenly spaced r-grid. The outer boundary 
condition is discussed later. 

The choice of an evenly spaced r-grid results in an unevenly spaced IP- grid. In essence, 
the IP-grid is stretched, but in a manner that is natural to the !P-grid, with more grid points 
at the outer edge where a finer grid spacing is desired. From equation (13), 


dlP = rpu dr. 


(18) 


Integration of equation (18) across one grid spacing we obtain 

AV k = «P* - = \[(rpu) k + (rpu)k -\]Ar (19) 

after using the trapezoidal rule. One consequence of equation (19) is that u cannot be zero 
outside the jet or d*P would go to zero. To get close to quiescent outer conditions, the 
free stream velocity is set to 1% of the initial highest velocity. 
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Figure 4. The stream function grid for mean flow calculations. 


Writing equation (14) as 


du d_ du_ 

dx ~ dY l A dW J’ 


( 20 ) 


where 


A = r 2 pufi e ff, (21) 

we implicitly difference this momentum equation, using the grid shown in Figure 4, as 
follows: 


U/+ u ~ UjM 

Ax 


AV k + x + AV k 


Ak + 


Uk + \ — Uk 


1/2 


AW k + 


Ak- 


k- 1/2 


U k ~ U k - I 


o+ 1 


( 22 ) 


The coefficient A is the average of adjacent grid values 

Ak + \n — \(Ak+Ak + i)' (23) 

Collecting terms, equation (22) may be arranged into tridiagonal form 

{ — Ak Ak - 1 / 2 w/r - 1 4* ( 1 + Ak Ak + 1/2 + A^ Ak - \a)iik — At Ak + 1,2 w* + 1 }; + 1 — «/,* » (24) 

where 


Ar = 


2Ax 

(AW k + i + AW k )A^r 


At = 


2 Ax 

(A'Fj' + t + A'F k )A'F k 


(25, 26) 
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Equation (15) is differenced in the same manner. The resulting tridiagonal equation is 
{ — A k D k - 1 / 2 w* .. i — A k B k . i ,2 ft k - i + (Am D k t- 1.2 + A k D k - ( :)i7* + ( 1+ A k A t- 1 : 

+ A k B k ~ \a)H k — A k D k + \afi k + 1 — A k B k + 1 2 H k + 1 }/ + 1 — Hi.k, (27) 

where 


B = rpu D = rpu 2 \ p eff - 


Pr, 


Petr 

Pr,« 


(28, 29) 


2.4.1. Boundary conditions 

Since the problem is axisymmetric, a symmetry boundary condition is enforced at r = 0, 
k — 1; that is, 

du/dV = dH/dV = 0. (30) 

The right side of equation (20) is differenced using half the grid spacing at the edge (see 
Figure 4): 


Uj+ 1,1 ~ Uj,\ 1 

Ax ~ AVzl2 



dii 

d¥ 


+ 1.1 + 1/2 



(31) 


From the boundary condition (30), the (j + 1, 1) term on the right side is zero. The final 
tridiagonal form is 


(Y, 2 Ax 2 Ax _ | 

+ ~aWI a ' + U2 J u ' ~ AWj A ' + t/2U2 j ~ UjA ' 


Similarly, for the energy equation we obtain 


(32) 


2 Ax 


2 Ax 


D 1 + \ t iU\ + ( 1 + £ y /2 + 1/2 )H 1 — D\ + 1/2M2 — 


2 Ax 


2 Ax 




B , 


AVI' 


+ 1/2 a ! — 
j } + 1 


(33) 


The outer boundary condition is simply that the u and H values equal the free stream 
conditions; that is, 

My+l./V+l — Uj.N + 1 — fix, Hj+ I.Af+l = ffj.N+t — Hx. (34) 

This also may be viewed as saying that derivatives in the jx-direction are zero at the outer 
edge. The consequence of this type of boundary condition is that the flow, expanding 
outward due to diffusion, must never reach the outer grid boundary or the calculations 
will be in error. The method by which this difficulty is overcome is discussed in the next 
section. 


2.4.2. Numerical method of solution 

Equations (24), (27), (32), (33) and (34), create a 2 x 2 block tridiagonal system for k - 1 
to N + l, where N is the number of grid spacings. The system can be solved by standard 
block tridiagonal routines [22]. The coefficient matrices of the block tridiagonal system are 
unknown at the j + 1 location since they are functions of the unknowns, fij+ u and 
Consequently, iterations must be made to complete the solution for each axial step. The 
procedure is as follows. 
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1. Guess u, f a and fl, • 11/ = 1, then use the initial conditions: 

Uz.k = M|.*, //:.* = //,*. 

For j > 1, we can extrapolate from previous solutions: 

m, + 1 .* = 2u,. k — M/ - 1 .* , Hj + 1 .* = 2//,. * — /// _ i.* . 

2. Solve for j + 1 values of p , m^/ and Pr,.,/-. 

3. Calculate d *£%!.* grid using equation (19). 

4. Solve the block tridiagonal system for new values of it J + t . k and H, + t .*. Also, new 
values of p, p eff and Pr <Y , are calculated. 

5. If the difference between the new solution and the previous solution is smaller than 
some convergence criterion, then the axial step is completed. Otherwise, the iteration 
process continues at step 3. 

To ensure, as mentioned previously, that the outer boundary condition is not 
compromised by the expanding flow, the grid must be large enough in the radial direction 
to encompass the flow. Rather than making the grid so large that the flow is encompassed 
at all axial locations, the outer edge of the initial grid is much closer to the outer edge 
of the jet, typically the grid ends at twice the outer radius of the initial jet. As the flow 
expands downstream, the outer edge of the jet flow is tracked until it is within some 
arbitrary distance from the outer edge of the grid; at which time, more grid points are 
added to the solution. All of the variables at the new grid points take on the free stream 
values that are assumed to always exist outside the grid. As the shear layer of the jet 
expands, it becomes possible to increase the dr-grid spacing since less grid points are 
necessary to define the shear layer. Thus, we do not simply continue to add grid points 
as the flow expands but, from time to time, we reduce the number of grid points by 
increasing the grid spacing. Actually, it is best to simply double the grid spacing so that 
the extrapolated guesses for the new variables occur along constant grid lines in the 
;t-direction. Thus, no interpolation is necessary. 

2.4.3. Numerical accuracy 

The numerical scheme for the governing parabolic equations, equation (22) for example, 
is first order accurate in both Ax and A W. The latter accuracy is due to the uneven spacing 
in the IP-grid. Since the mean flow results are used in the stability calculations, it is found 
there that to obtain converged stability results a finely spaced r- grid is necessary for the 
calculation of the derivatives of the mean flow variables. Using a fine r- grid results, 
through equation (19), in a fine, but unevenly spaced V'-grid. An example to illustrate the 
convergence of the mean flow calculations using the mixing length turbulence model is 
shown in Figure 5 for three different grid spacings. The calculated shear layer half-width 
b and the vorticity thickness <5 are shown for the shear layer of a single jet under both 
incompressible and compressible conditions. The half-width shows the ability of the code 
to find the defined edges of the shear layer in order to calculate the half-width and vorticity 
thickness shows its ability to calculate accurately the maximum derivative of the mean 
velocity profile. For the range of grid spacings shown, the mean flow results are 
independent of the grid spacing. 


3. NUMERICAL PREDICTION 
3.1. SINGLE JET CALCULATIONS 

The numerical predictions are compared in Figure 6 to data taken form Zaman [33] for 
a subsonic, M — 0-5 axisymmetric jet. The axial variation of the centerline velocity is 
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Incompressible single jet parameters; (b) compressible single jet parameters. 


shown in Figure 6(a). The calculations show good agreement with the data and to an 
equation for centerline velocity from Witze [34], based on a correlation of many different 
jet measurements. In Figure 6(b) are shown the growth of the jet half-velocity width and 
the ratio of the jet half-velocity width to momentum thickness. The calculations 
underpredict the spread of the jet as defined by the half-velocity width. 

Results for higher speed jets are shown in Figure 7. For the M = 1-37 jet, the centerline 
velocity data are taken from Lau et al. [6] (Figure 7(a)) and the half-velocity point data 



Figure 6. A comparison of subsonic single jet calculations (Af = 0-5) to measured data [33]. (a) Centerline 

velocity: , calculated; O, data; , Witze [34] centerline; (b) jet half-velocity width/jet exit diameter, □ 

(left scale); jet half-velocity width/momentum thickness. A (right scale); , calculated. Operating conditions: 

t/, = 172-5 m/s, T, = 295-4 K. Ri = 1-27 cm. 
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Figure 7. A comparison of subsonic single jet calculations to measured data, (a) centerline velocity: (b) jet 

half-velocity point. Operating conditions: M = 2-22, U\ - 538 m/s, T x = 147 K, R\ = 1-28 cm; , calculated; 

□ , data [36]; M = 1-37, U\ = 472 m/s, T\ = 295-4 K, R\ = 2-55 cm; , calculated; A. data [6, 35]. 


is from Lau [35] (Figure 7(b)). The calculated core region spreads at the measured rate, 
although the measured core length is slightly shorter than predicted. For the M = 2-22 jet, 
the comparisons for both the centerline velocity and the jet half-velocity point are made 
with data from Eggers [36]. Again, the calculations underpredict the measured spreading 
of the jet. These results show the inaccuracies that exist in trying to model the effects of 
compressibility in our mean flow calculations. The calculated spreading rate of a 
compressible jet is determined by the compressible calibration equation (A3) given in 
Appendix A. The measured spreading rate data are scattered around this correlation and 
the predictions reflect these discrepancies. Any measured jet spreading rate data that are 
larger than that predicted by the calibration equation will have calculated spreading rates 
that are less than the measured spreading rates. This is reflected in the comparison between 
the calculated jet half-velocity point and the data measured by Eggers. Since the 
experimental jet spreading rate is higher than that predicted by the calibration equation, 
the measured jet half-velocity point moves outward faster downstream than predicted 
(Figure 7(b)). Experimental jets that have spreading rates much closer to the calibration 
curve (e.g., the measured Lau data) have better agreement between the calculated and the 
measured jet half-velocity points. As with many other attempts at turbulence modelling, 
the results are only as good as the calibration and it is difficult to apply the same model 
to experiments for a wide range of operating conditions in many different facilities. Also, 
the turbulence model is calibrated to match growth rates in the initial core region of a jet, 
the region of most importance to growing instabilities. In the downstream region, the jet 
is a thick axisymmetric shear layer and the current turbulence model, based on thin shear 
layer results, provides less representative growth rates. 

A final comparison is made for a single jet in a moving stream. These results are 
important to the subsequent consideration of coaxial jets on the effect of Mach number 
and the effect of velocity ratio on jet spreading characteristics. The effect of the velocity 
ratio on the centerline velocity is shown in Figure 8(a). The calculated results show that 
the turbulence model gives a decrease in jet spreading and an increase in the potential core 
length as the velocty ratio increases from 0 096 to 0-497. This agrees with the results 
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measured by Morris [29] for the effects of velocity ratio on centerline velocity. The Mach 
number effect at a velocity ratio of about 0- 1 is shown in Figure 8(b). Again, the calculated 
results behave in the same manner as the measured data as the Mach number increases 
from 0-47 to 1-67. Hence, the turbulence model allows for the calculated jet behavior to 
mimic the expected jet behavior for a compressible jet in a co-flowing stream. 

3.2. COAXIAL JET CALCULATIONS 

With the single shear layer turbulence model tested on both subsonic and supersonic 
jets with and without a co-flowing stream, we now add an outer jet flow to create a coaxial 
jet with two shear layers. Both NVP and IVP jet examples are considered. A comparison 
of calculated velocity profiles for a NVP jet to measured data [30] at four axial locations 
is shown in Figure 9. Initially, the agreement is very good. The outer shear layer is 
spreading faster than the inner shear layer, as expected from the results in Figure 8 for 
a jet in a co-flowing stream. Downstream, the calculated jet centerline velocity begins to 
decay sooner than the measured data. This again reflects the differences between the 
correlation for the spreading rate and the measured values for this set of data. However, 
in essence, the calculations produce a reasonable representation of the merging jet for an 
initial normal velocity profile. 

Next, the inverted velocity profile data from von Glahn et al. [37] are used for 
comparison with calculations. Velocity profiles measured at five axial locations are shown 
in Figure 10. The IVP predictions show similar velocity profiles as the two shear layers 
spread and merge downstream. The other calculated results shown in the figure are from 
Georgiadis and Yoder [20]. Their results are indicative of the problems encountered when 
using Navier-Stokes solvers with k-e turbulence models. The initial profiles are not smooth 
and the spreading of the shear layer downstream is greatly underpredicted. 

A final comparison with IVP jet data is shown for both velocity and temperature profiles 
in Figures 11 and 12: the data are from Tanna et al. [38]. In this case, the outer stream 
is much thinner than in the previous example causing the two shear layers to mix sooner. 



Figure 8. A comparison of single jet in moving stream calculations to measured data [29]. (a) Centerline 

velocity variations with external flow: , calculated, M =0-47, VijV i = 0 096; O. data; , calculated, 

M = 0-47, Ui/Ui =0-497; [>, data, (b) Centerline velocity variations with Mach number: , calculated, 

M = 0-47. t/’/t/i = 0-096; O. data; . calculated, M - 1 -67, = 0-098; O, data. Operating conditions 

R t = 2-54 cm; M = 0-47, £/, = 158 m/s, Ti = 283 K; M = 1-67, £/, = 461 m/s, T\ = 190 K. 
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Figure 9. A comparison of normal velocity profile jet calculations to measured data [30]: , calculated; O, 

data. Radial velocity profiles: (a) x = 4; (b) x = 8; (c) x - 16; (d) x = 32. (e) Centerline velocity. Operating 
conditions: U, = 411 m/s, T t = 657 K, R t = 1-96 cm; Ui = 274 m/s, T 2 = 292 K, R 2 = 3-91 cm. 





Figure 10. A comparison of inverted velocity profile jet calculations to measured data [37]: , calculated; 

O. data; , calculated in reference [20]. Radial velocity profiles: (a) x = 2; (b) x = 4; (c) x = 8; (d) x = 16. 

(e) x = 26. Operating conditions: U\ = 314 m/s, T t = 248 K, R t = 5 03 cm; U 2 - 459 m/s, T 2 = 1029 K, 
Ri = 8-80 cm. 











COAXIAL JET MEAN FLOW 


659 



0 1 2 3 401 2 3 4 


r 

Figure 11. A comparison of inverted velocity profile jet calculations to measured data [38]: , calculated; 

O, data. Radial velocity profiles: (a) x ~ 0-38; (b) x - 2-7; (c) x = 3-8; (d) x = 1 1 -5; (e) x — 23 0. Operating 
conditions: Ut = 273 m/s, T\ = 434 K, R\ = 2-61 cm; Ui = 477 m/s, Ti = 758 K, R 2 — 3-57 cm. 



r 


Figure 12. A comparison of inverted velocity profile jet temperature calculations to measured data [38]. See 
Figure 1 1 for legend and operating conditions. 
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Both good qualitative and quantitative agreement is seen for the velocity and temperature 
profiles. Thus, the present calculations are capable of producing a reasonable merging jet 
with an initial inverted velocity profile. 


4. SUMMARY 

Much of the thoretical work on the stability of free shear flows relies on the analytical 
definition of the mean flow in order to solve for the stability characteristics. Experimental 
measurements of shear flows have provided the basis for simple analytical forms to 
describe the mean flow in the stability calculations. Such mean flow data are lacking for 
two stream coaxial jets, especially at supersonic conditions. The data that are available 
are not sufficient to develop analytical expressions for the mean flow properties of coaxial 
jets that transition smoothly from two streams to one fully developed jet downstream for 
both normal velocity profile and inverted velocity profile cases. A numerical scheme has 
been developed to calculate the mean flow properties of a coaxial jet with two initial shear 
layers. The development of a turbulent, axisymmetric free jet is assumed to be governed 
by the Reynolds-averaged, compressible, parabolic boundary layer equations. To complete 
this set of equations, a turbulence model is necessary. Any turbulence model chosen for 
these calculations would have required modification and calibration to obtain reasonable 
results for both normal and inverted profile coaxial jets. For expediency, a modified 
and calibrated mixing length model has been used to represent the turbulent stresses 
and obtain an effective viscosity. Mean flow properties have been predicted for coaxial 
jets with either a normal initial velocity profile or an inverted initial velocity profile 
that expand and merge in a manner representative of measured data. Both good qualitative 
and quantitative agreement has been found for both velocity and temperature profiles. 
These mean flow results can be used to calculate instability wave characteristics for coaxial 
jets. These calculations are described in subsequent papers for both NVP [4] and IVP [5] 
jets. 
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APPENDIX A: CALIBRATION OF THE MIXING LENGTH MODEL 


A. I. INCOMPRESSIBLE FACTOR Ci 

For given flow conditions, the mixing length model must provide the effective viscosity 
that causes the jet shear layer to spread. Ci is the only factor available for adjustment in 
order to obtain the proper spreading rate when the flow is incompressible. Many authors 
have used experimental evidence and arguments about large scale structures dominating 
the development of turbulent shear layers (some recent references are references [28] and 
[39-43]) to develop a prediction for the incompressible vorticity thickness growth rate S,Jx 
as a function of the velocity ratio r and the density ratio s. 


— = 0-088 + V S \ 

x 1 +rjs 


(Al) 


where r = U 2 /U\ and s — pi/pi- The constant 0-088 used here is an average of the constants 
from the above references when all are converted to a vorticity thickness growth rate, either 
from information given in the reference or from the assumption of an error function shape 
for the velocity profile. Although equation (Al) has been developed for a planar shear 
layer, it is assumed to hold in the initial core region of a jet where there is a constant 
velocity at both edges of the shear layer. 

The calibration of Ci is carried out for a range of expected r and s values. Given r and 
s, Ci is varied until the calculated vorticity thickness growth rate in the core region agrees 
with that predicted by equation (Al). To ensure that p T dominates the effective viscosity 
and that the flow is essentially incompressible, the calculations have been carried out on 
a low speed, large diameter jet such that the Reynolds number was greater than 10 6 and 
p T »p. The resulting family of curves is correlated to develop an equation for C t for r < 1 
and an equation for r > 1. The curves take the form 


where, for r < 1, 


and, for r > 1, 


Ci(r,s) 


A(r) + B(r)y/s -F C{r)s 
(1 + ^)(1 + r yj~s) 


A{r) = (6-5919 + 11-918/- - 4-1855/- 2 ) x 10“ 2 , 
B(r) = (10-880 - 2-3578r + 6-5642r 2 ) x 10~ 2 , 
C(r) = (3-1013 4- 16-420r - 5-8217r 2 ) x 10~ 2 


A(r) = (101-95 + 51-507/- - l-6000r) x 10~\ 
B(r) = (60-287 + 75-160/- + 1-3818H) x 10“ 3 , 
C(r) = (723-97 + 753-34r - 6-2101/- 2 ) x 10~ 4 . 


(A2) 


A. 2. COMPRESSIBILITY FACTOR C 2 

To develop the compressibility factor, experimental results for compressible shear layers 
have been used to make a plot of normalized vorticity thickness growth rates versus 
convective Mach number [12]. A simple curve fit to this data is given by 


S;„(M c )/SL = 1 + 0-785[exp ( — 2 A/ 2 ) - 1], 


(A3) 
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where the convective Mach number is 


M r = M\ 


S*\ 



(A4) 


S' m , is the incompressible spreading rate for the same r and s , and A/, is the Mach number 
for the stream with velocity U } . 

For given r, s and M\, C 2 is adjusted, allowing C\ to take on its value defined by equation 
(A2), until the calculated compressible vorticity thickness growth rate in the core region 
agrees with the value predicted by equation (A3), when d' 00 is given by equation (Al). The 
results for C 2 have been correlated with M c to obtain the following equation: 


C 2 (M e ) = 1 + 04959[exp (-1*4593 M; + 0-0427 M] - 0-3658M t 4 )- 1], (A5) 



